Validity of the Wigner-Seitz approximation in neutron star crust 
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Since the seminal work of Negele and Vautherin, the Wigner-Seitz approximation has been widely 
applied to study the inner crust of neutron stars formed of nuclear clusters immersed in a neutron 
sea. In this article, the validity of this approximation is discussed in the framework of the band 
theory of solids. For a typical cell of 200 Zr, present in the external layers of the inner crust, it is shown 
that the ground state properties of the neutron gas are rather well reproduced by the Wigner-Seitz 
approximation, while its dynamical properties depend on the energy scale of the process of interest or 
on the temperature. It is concluded that the Wigner-Seitz approximation is well suited for describing 
the inner crust of young neutron stars and the collapsing core of massive stars during supernovae 
explosions. However the band theory is required for low temperature transport properties as, for 
instance, the effective neutron mass giving rise to entrainment effects. 
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In the standard model of neutron stars the crust 
is believed to be formed of nuclear clusters in a body 
centered cubic lattice stabilized by Coulomb forces and 
considered infinite and pure (made of only one type of 
nuclei at a given density). In the inner crust, at densi- 
ties above ~ 4.10 11 g.cm -3 and below ~ 10 14 g.cm" 3 , the 
"neutron drip" regime is reached and the clusters are sur- 
rounded by a neutron fluid. A formal comparison can be 
made with electrons in ordinary solids present on earth: 
part of the neutrons participate to the nuclear clusters 
which form the lattice (equivalent to electrons bounded 
to atoms) while part of the neutrons are delocalized over 
the whole crystal (equivalent to valence electrons). As a 
consequence, the band theory of solids developed in con- 
densed matter [2[ can be applied to describe the crust 
of neutron star. But due to the highly specific numer- 
ical issues of band theory, nuclear physicists have pre- 
ferred to use an approximation due to Wigner and Seitz 
(W-S) [1, 0|, where the crust is divided into indepen- 
dent and spherical cells. Since the work of Ncgclc and 
Vautherin |5|, the W-S approximation has been used to 
predict the structure of the crust, the pairing properties, 
the thermal effects, or the low lying energy excitation 
spectrum @, 0, H ©, EE [U G3 . Only recently, band the- 
ory calculations have been carried out in order to study 
the hydrodynamical properties of the neutron fluid and 
in particular the neutron effective mass giving rise to en- 
trainment effects [U, [3 EH , although these calculations 
are not yet self-consistent. While the W-S approxima- 
tion is well justified below the "neutron drip" regime, its 
validity beyond remains to be assessed. 

In this article, we investigate the limitations of the W- 
S approximation in the p ~ 7.10 11 g.cm -3 density layer 
of the inner crust, composed of a crystal of zirconium like 
clusters [|| surrounded by the neutron gas. In Sect.Q] be- 



fore discussing the W-S approximation, we briefly review 
the band theory of solids. Then we compare in Sect. |TT] 
the results of the band theory with those of the W-S 
approximation for the single particle wavefunctions and 
energy spectra. Consequences for the properties of the 
neutron gas are discussed. 



I. MICROSCOPIC QUANTUM DESCRIPTION 
OF NEUTRON STAR INNER CRUST 

An accurate description of the inner crust, assuming 
that it is a perfect crystal, should rely on the band theory 
of solids [|J . In this section, we briefly review this theory 
in the context of neutron star crust and discuss the W-S 
approximation in this framework. 



A. Band theory of solids 

According to the Floquet-Bloch theorem, the single 
particle quantum states are given by modulated plane 
waves 



<Pak(r) = u ak (r)e l 



k r 
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where the functions u a k(r) have the full periodicity of 
the lattice. Each single particle quantum state is thus 
labeled by a discrete index a and by a wave vector k. 
The energy spectrum is therefore formed of a series of 
sheets or "bands" in fc-space. 

The Bloch states ([T]) are completely determined by the 
knowledge of the functions u a k(r) inside a W-S cell of 
the lattice, whose shape is imposed by the symmetry of 
the crystal. The cell, centered around one nuclear clus- 
ter, is electrically neutral and therefore contains as many 
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electrons as protons. The effects of the ion lattice on 
the electrons, which give rise to complicated band struc- 
tures in ordinary terrestrial matter are negligible in the 
inner crust of a neutron star due to the very high den- 
sities [161] . Nevertheless, the neutron band effects due to 
nuclear inhomogeneities cannot be ignored. 

In the present study, we consider the outermost layers 
of the inner crust where pairing effects are negligible [17| . 
In the Hartree-Fock approximation with Skyrme forces 
which we shall consider in the following, the occupied 
nucleon single particle wave functions are obtained by 
solving the self-consistent equations (q = n,p for neu- 
trons and protons respectively) 

h { o q) &r)=e%^(r) (2) 
where the single particle Hamiltonian is defined by 



2mf(r) 



V + U q {r) -iW q {r) - V xa, (3) 



the effective masses m®(r), mean fields U q (r) and spin- 
orbit terms W q (r) being functional of the single particle 
wave functions. These equations have to be solved inside 
the W-S cell with the boundary conditions imposed by 
the Floquet-Bloch theorem 



(4) 



where T is any lattice vector. This means in particular 
that the wavefunction between two opposite faces of the 
cell has a phase shift e where T is the corresponding 
lattice vector. The single particle energies are periodic in 
the reciprocal lattice whose vectors K satisfy K T = 2im 
(where n is any integer) 
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Consequently only the values of k inside the first Brillouin 
zone (i.e. W-S cell of the reciprocal lattice) are relevant. 
Equivalcntly, equations (T2J) can be written directly for 

the u^l (r) functions in the decomposition (p} which leads 
to 
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where the fe-dependent Hamiltonian Yi^f 1 is defined by 



h 2 k 2 
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hk . 



(7) 



and the velocity operator v q is defined by the commutator 



= — \r h [ ' 



(«)i 



(8) 



The band theory takes into account all the symmetries 
of the system. However equations ^ with the boundary 
conditions Q are numerically very complicated to solve. 
The approximation introduced a long time ago by Wigner 
and Seitz in the study of metallic sodium [J H[ has been 
widely applied in the context of neutron star crust, as 
described below. 



B. Wigner-Seitz approximation 

The spherical W-S approximation is a computation- 
ally very efficient method with the advantage of reduc- 
ing the 3D partial differential Eqs. ([2]) to ordinary differ- 
ential radial equations. This approximation is twofold. 
First of all, the Hamiltonian /il in equation ^ is ne- 
glected. Consequently the wave functions and the en- 
ergies are independent of k and approximated by the 
solutions at k = 0. Only the band index a remains. 
Secondly, the W-S polyhedron is replaced by a sphere 
of equal volume. The equations are then usually solved 
with the Dirichlet- Neumann mixed boundary conditions 
which yield a nearly constant neutron density outside the 
cluster. 

The W-S approximation turns out to be very good if 
the boundary conditions play a minor role. For instance, 
bound states whose associated wave functions are van- 
ishingly small outside the clusters are very well treated 
provided that the spatial extent of these states is smaller 
than the lattice spacing. This condition is fulfilled al- 
most everywhere in the crust except in the bottom layers 
where the clusters nearly touch. The aim of this paper is 
to investigate the validity of the W-S approximation for 
the outermost layers of the inner crust where the bound 
neutron states are not altered by the boundary condi- 
tions. 

Let us emphasize that in the W-S approximation, the 
nuclear clusters are supposed to be spherical while in 
the full band theory, no assumption is made about their 
shape. For the low densities of interest in this study, the 
nuclear clusters can still be considered as spherical. It 
should mentioned that in a recent development of the 
W-S approximation (lgj |. the W-S cell is replaced by a 
cube with strictly periodic boundary conditions. Possi- 
ble deformations of the nuclear clusters are thus included 
but at the price of unphysical boundary conditions be- 
cause the W-S cell of the body centered cubic lattice is a 
truncated octahedron and not a cube (the cube being the 
W-S cell of the simple cubic lattice). This is the reason 
why we still consider the spherical W-S approximation 
closer to the physical situation than the cubic one at low 
density. 



II. COMPARISON BETWEEN THE BAND 
THEORY AND THE W-S APPROXIMATION 

The comparison between the band theory and the W-S 
approximation gives an estimate of the contribution ([7]) 
of the fc-dependent Hamiltonian hk and incidentally on 
the effects of the boundary conditions. We have consid- 
ered the shallow layers of the crust at the average baryon 
density p ~ 7 x 10 11 g/cm 3 formed of a crystal made 
of zirconium like clusters (Z — 40) with 160 neutrons 
(bound and unbound) per lattice site||. In the follow- 
ing we shall refer to such cluster as 200 Zr. Under such 
conditions, only unbound neutrons are sensitive to the 
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FIG. 1: Neutron (full lines) and proton (dotted lines) density 
distributions for 200 Zr in the W-S cell, obtained in the W-S 
approximation. 



FIG. 2: Unbound neutron density calculated with the W-S 
approximation (WS, full lines) and the full band theory (BT, 
dotted lines). 



boundary conditions and the nuclear clusters are spheri- 
cal. 

For the comparison, we first solve the self-consistent 
Hartree-Fock equations in coordinate space, considering 
the W-S approximation. The effective force used is the 
same as in references @,|8, 9], namely the Skyrme interac- 
tion SLy4 fl9j | . In Fig.[TJ we show the neutron and proton 
densities calculated in the spherical W-S cell with 200 Zr. 
The size of the box is i? C cii = 49.2 fm. The cell exhibits 
a very large and diffuse neutron skin, which is typical of 
those systems pj. A small but non-zero neutron density 
is present at large radius generating a non-zero mean field 
potential. Asymptotically, this potential is equal to -0.05 
MeV. All the states with energy larger than -0.05 MeV 
are therefore unbound or "free". We found that among 
the 160 neutrons per lattice site, 70 are unbound. 

The effective mass m®(r) and the mean field potential 
U n (r) obtained for the spherical cell are used to con- 
struct an effective Schrodinger equation for band theory 
calculations. As the spin-orbit splitting is weak for most 
of the states (see Fig. [SJ, we set the spin-orbit poten- 
tial W n (r) to zero. In order to study the effects of the 
boundary conditions, the Schrodinger equation is solved 
with no further iterations by imposing the Bloch bound- 
ary conditions (U]) and u sing the Linearized Augmented 
Plane Wave method (see [ljf for details). The coordinate 
space is divided in two regions: a spherical region around 
the cluster plus an interstitial region. In the latter region, 
the wave functions are expanded on a plane wave basis in 
order to fulfill the Bloch boundary conditions. The lat- 



tice spacing is determined by requiring that the volume 
of the W-S sphere is equal to the volume of the exact 
W-S cell of the crystal, assumed to be a body centered 
cubic lattice [j]]. 

In the following, we compare the single particle wave 
functions and energy spectra of the unbound neutrons. 



A. Single particle wavefunctions 

As already discussed in Sect. II B[ the wave functions 
of the bound states are nearly independent of the bound- 
ary conditions. As a consequence, we expect that the 
band theory and the W-S approximation provide identi- 
cal bound states. Since the unbound states are orthog- 
onal to the bound states, the W-S approximation and 
the band theory are expected to yield similar unbound 
wavefunctions inside the nuclear clusters. This is con- 
firmed by the calculation of the density distribution of 
the unbound neutrons (whose single particle energies ex- 
ceed —0.05 MeV as discussed previously) shown on Fig. [21 
obtained with the band theory and the W-S approxima- 
tion. For the comparison, the density p(r) obtained from 
the band theory, has been averaged over the solid angle 
around one lattice site as follows 

p(r) = J^p( r )> ( 9 ) 

where r is the radial distance from the lattice site. Simi- 
lar density oscillations are obtained in both calculations 
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in the vicinity of the nuclear cluster for r < 10 fm as 
expected. Qualitative differences in the neutron density 
are however observed in the interstitial regions outside 
clusters due to different boundary conditions. The un- 
bound neutron density distribution is nearly flat in the 
band theory while it is more fluctuating in the W-S ap- 
proximation. The Bloch wavefunctions outside the clus- 
ters are similar to plane waves (thus giving a constant 
density) which cannot be properly described in the W-S 
approximation owing to the expansion of the wavefunc- 
tions into only a few spherical harmonics. An analysis of 
the contribution of each single particle wave function to 
the unbound neutron density in the W-S approximation 
reveals that the oscillations at small radius are mainly 
coming from p-states, such as 3pi/2 or 3p3/2 whereas at 
larger radii (r > 20 fm) only a few larger I states, mainly 
d-f-g-h states, are contributing to the free neutron den- 
sity. 

As a result, the W-S approximation predicts a different 
number of neutrons outside the cluster than in the band 
theory. Since the total number of free neutrons should 
be the same in both calculations, the difference in the 
density profile at large radius implies a larger difference 
in magnitude at small radius. This is more clearly seen 
on Fig.[3j by plotting the integrated number N(r) of free 
neutrons inside the W-S cell at radius r, defined by 

N(r) = 4tt / r' 2 p(r')dr' , (10) 
Jo 

p(r) being the local density of unbound neutrons. The 
figure shows that the W-S approximation underestimates 
the number iVj n = N(R) of free neutrons inside the re- 
gion of radius R around the cluster and consequently 
overestimates the number 7V out = N(R CC \\) — N(R) of 
free neutrons outside. For both calculations, the number 
of free neutrons in the cell is N(R ce a) — 70. Quantita- 
tively taking R = 15 fm, the difference between the two 
calculations is about AN = |AJVjj,| = |AJV ou t| = 3 which 
is rather small. 

This first comparison shows that the single particle 
wave functions of unbound neutrons are qualitatively well 
reproduced by the W-S approximation inside the nuclear 
clusters. The main differences in the wavefunctions be- 
tween the two calculations are found in the interstitial 
region due to the different boundary conditions. How- 
ever this has a rather small effect on the ground state 
properties of the neutron gas, like the neutron density 
distribution. More generally the W-S approximation can 
be expected to be a good approximation to the full band 
theory for evaluating the matrix elements of any operator 
taking vanishing values outside the cluster region. 

B. Single particle energy spectrum 

Figs. 2] and O show the energy spectrum of the unbound 
neutrons obtained in the band theory and in the W-S ap- 
proximation, respectively. It should be noted from Fig. [5] 
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FIG. 3: Integrated unbound neutron number (see text) cal- 
culated with the W-S approximation (WS, full lines) and the 
band theory (BT, dotted lines) 




FIG. 4: Unbound single particle energy spectrum as obtained 
in the band theory vs the Bloch wave vector k, along high 
symmetry lines in the first Brillouin zone using standard no- 
tations 20] . 



that the spin-orbit splitting is very weak for d, f, g and 
h states (as predicted by Negele&Vautherin ||) but not 
for p states. This is due to the fact that the spin-orbit 
splitting is proportional to the convolution of the density 
gradient together with the wave functions. The density 
gradient is localized in the central cluster while the d, f, g 
and h states are mostly in the external region. For those 
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FIG. 5: Unbound single particle energy spectrum obtained in 
the W-S approximation. 
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states, the convolution leads to a weak splitting. 

Since in band theory, the energies depend also on the 
wavevector fc, only the energy bands along some specific 
symmetry directions in A;-space are displayed. The en- 
ergy spectrum obtained with the W-S approximation is 
comparable to the one obtained in the band theory for 
the symmetry point T corresponding to the center k = 
of the first Brillouin zone (W-S cell of the reciprocal lat- 
tice) . The correspondence is not exact and the differences 
come from the spherical approximation. As a result, the 
W-S approximation predicts less states but with larger 
degeneracies than the band theory at the symmetry point 
r. The figures show clearly another important difference 
between the band theory and the W-S approximation: in 
the former case the energy spectrum is continuous while 
in the latter case it is discrete. 

A relevant quantity to compare the energy spectra is 
the level density, which plays a pivotal role when calcu- 
lating dynamical processes, ft is defined by 



9(E) 



Z Jbz 



d 3 fc 



5{E-e ak ) 



(11) 



where the integral is taken over the first Brillouin zone 
(BZ). ^ 

Using the 8 function to integrate out one of the vari- 
ables, the level density becomes 



(2tt) 3 



(12) 



where the integral is taken over the surface of constant 
energy e a f. = E in fc-space. Expression (TT21 shows that 
the level density is a probe of the topology of constant en- 
ergy surfaces in A;-space. We have extracted the level den- 
sity from band theory by using the Gilat-Raubenheimer 
method as in reference flil [l5j . 



FIG. 6: Level density for neutron unbound states calculated 
with the band theory (solid line) compared with the predic- 
tion of the free Fermi gas (dashed line). The energy resolution 
is of the order of keV. 



In the W-S approximation the level density reduces to 
a discrete sum 

g ws (E) = J2( 2 j + 1 )S(E-£ nj e). (13) 

njl 

In Fig. [6l we show the level density predicted by the 
band theory for the unbound single particle levels. As ex- 
pected the figure shows that the energy spectrum in the 
band theory is continuous and has a complex fine struc- 
ture. The spectrum exhibits a quasi band gap of about 
30 keV slightly around 90 keV. This is in sharp contrast 
with the W-S approximation for which the energy levels 
are discrete and separated by about 100 keV. In other 
words, the W-S approximation overestimates the neu- 
tron shell effects. The global energy dependence of the 
level density follows the behavior of the free Fermi gas, 



1(E) 



Vcell 
2^" 



E v , 



(14) 



where E v ~ —0.031 MeV is the energy at the bottom 
of the valence band. The agreement between the two 
curves is very good for energies close to E v . This means 
that the Fermi surface is nearly spherical at low ener- 
gies as confirmed by the calculations shown on Fig. [7] 
This is due to the fact that the Fermi wavelength of the 
unbound neutrons is much larger at low energies than 
the lattice spacing. As a consequence the effect of Bragg 
diffraction is negligible. It can be inferred from Fig. [5] 
that distortions of the Fermi surface from the spherical 




FIG. 7: (Color online) Constant neutron energy surface of a 
body centered cubic lattice of 200 Zr inside the first Brillouin 
zone for E = -0.0015MeV. 

shape happens at energies larger than 0. The first kink 
around the zero of energy is a characteristic van Hove 
singularity (as a result of the vanishing of the gradient 
Vfc s a k at some k points in the expression (|12p. also vis- 
ible on Fig. 2]) and indicates a topological transition of 
the Fermi surface. This occurs when the Fermi sphere 
touches the faces of the first Brillouin zone. For a body 
centered cubic lattice, the transition takes place when 
the radius of the sphere is equal to y/2w/a where a is the 
lattice spacing. Above the first kink, the Fermi surface 
becomes non spherical with the appearance of necks close 
to the Brillouin zone faces as illustrated on Fig. [H] The 
Fermi surface undergoes further topological changes as it 
crosses Bragg planes (higher Brillouin zones) as revealed 
by the singularities in the level density. The actual Fermi 
surface (associated with the 160 neutrons per cell) has 
a very complicated shape with 11 branches (associated 
with the 11 bands which cross the Fermi level). 

In Fig. [9l we show the integrated number of single 
particle levels which are below a given energy E 

N(E) = [ deg(e). (15) 

The W-S approximation exhibits a stair steps structure 
due to the discretization of the states. As already no- 
ticed, the energy levels are discrete and highly degen- 
erate, due to the imposed spherical symmetry around 
the cluster. In contrast, in the band theory the spher- 
ical symmetry is partly broken due to the translational 
symmetry of the crystal lattice. The energy levels thus 
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FIG. 8: (Color online) Constant neutron energy surface of a 
body centered cubic lattice of 200 Zr inside the first Brillouin 
zone for E = 0.0085MeV. 



broaden into bands, with a low residual degeneracy (at 
most 6 at each k point for cubic crystals, counting the 
spin degeneracy 20]), which overlap so that the energy 
spectrum is continuous. 

In conclusion, for processes which involve transfered 
energies above the characteristic level spacing around 
the Fermi energy, the differences between the W-S ap- 
proximation and the full band theory are expected to be 
small. For instance, it is typically the case for neutrino 
response function [23] or thermal effects before the star 
cools completely down. However at lower energies, as 
pertaining for instance the effective neutron mass rele- 
vant for fluid dynamics [H, [3 [l5| , the full band theory 
is required. The level spacing can be roughly evaluated 
from the quantity h 2 /2mR 2 cU . From the top to the bot- 
tom layers of the crust, the characteristic level spacing 
varies from about 100 keV to 200 keV, which corresponds 
to temperatures of the order of 10 9 K. Such temperatures 
are found in young neutron stars less than a few hundred 
years after their birth [2^]. The W-S approximation is 
therefore well suited for describing the hot dense matter 
(except for the high density layers where the spherical 
approximation may be too restrictive) in young neutron 
stars and in the collapsing core of massive stars during 
supernova explosions [6j. This discussion however does 
not take into account pairing effects, which are negligible 
in the outermost layers of the crust considered in this 
work but are expected to be important in denser and 
deeper layers [13] • It should be noted that a recent study 
has shown that the pairing properties of the unbound 
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FIG. 9: Integrated state number calculated with the W-S 
approximation (WS, dashed lines), band theory (BT, solid 
line), and free Fermi gas (doted lines) 

neutrons are strongly sensitive to the choice of boundary 
conditions in the W-S approximation, especially in the 
bottom layers of the crust [12, HH ■ 

III. CONCLUSION 

In this article, a comparison has been done between 
the full band theory and the W-S approximation which 
has been widely applied in studies of neutron star crust. 
The external layers of the inner crust at a baryon density 
p ~ 7 x 10 11 g/cm 3 composed of zirconium like clusters 
200 Zr have been considered. Since the bound nucleons are 
not much affected by the boundary conditions, we have 
focused on the unbound neutrons. We have shown that 
the ground state properties such as the unbound neutron 
density distribution, are rather well reproduced by the 
W-S approximation, while the dynamical properties de- 



pend on the process of interest or the energy exchanged. 
It should also be noticed that depending on the quantities 
of interest, the free neutron model for the unbound neu- 
trons could be a good first approximation. In the future, 
it could be interesting to explore an intermediate scheme 
which goes beyond the Wigner-Seitz approximation and 
remains simpler to implement in numerical calculations 
than the full band theory. 

The energy spectrum is continuous in the full band the- 
ory with no energy gaps while the W-S approximation 
yields a discrete spectrum thereby overestimating neu- 
tron shell effects. The W-S approximation can therefore 
be applied whenever the processes under considerations 
involve energies larger than the level spacing induced by 
the discretization, which in the present case is of order 
~ 100 keV. In particular, the W-S approximation is well 
suited for describing the hot (T > 10 9 K) dense matter in 
the inner crust of young neutron stars and in the collaps- 
ing core of massive stars during supernovae explosions. 
However low temperature transport processes such as the 
effective neutron mass relating the momentum to the ve- 
locity and giving rise to entrainmcnt effects, require a 
fine knowledge of the energy spectrum around the Fermi 
level (i.e. the Fermi surface) which cannot be reproduced 
by the W-S approximation. Since the lattice spacing is 
predicted to decrease with increasing depth, becoming 
comparable to the size of the nuclear clusters and to 
the Fermi wave length of the free neutrons at the bot- 
tom of the crust, the validity of the W-S approximation 
should be carefully investigated in the denser layers of the 
crust, especially concerning pairing effects. Besides, the 
assumption of spherical symmetry in the W-S approx- 
imation is probably too restrictive near the crust-core 
interface where the clusters are expected to be strongly 
deformed [l6l |. 
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